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ABSTRACT 

This  document  describes  the  results  obtained  using  two  methods  for  ab  initio 
calculation  of  intermolecular  potential  parameters  for  gaseous  decomposition  products 
of  energetic  materials:  a  Multipole  Expansion  method,  suitable  for  axially  symmetric 
molecules,  and  a  Monte  Carlo  method,  which  can  be  used  to  obtain  temperature 
dependent  average  potential  energy  parameters  for  any  molecule. 
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Ab  Initio  Calculation  of  Intermolecular  Potential 
Parameters  for  Gaseous  Decomposition 
Products  of  Energetic  Materials 

Executive  Summary 

The  weapons  systems  of  the  Australian  Defence  Force  rely  on  energetic  materials  for 
rocket  and  projectile  propulsion  (propellants)  and  terminal  effects  (explosives).  In 
addition  to  research  aimed  at  decreasing  energetic  material  sensitivity  for  greater 
safety,  improvements  in  performance  are  also  being  sought  for  increased  range  and 
thrust  control  in  propellants  and  for  enhanced  terminal  effects  in  explosives. 

It  is  necessary  to  model  the  decomposition  of  energetic  materials  (combustion  for 
propellants  and  detonation  for  explosives)  in  order  to  gain  a  predictive  capability  for 
the  performance  of  existing  and  potential  energetic  materials  which  may  be  utilised  in 
future  weapons  systems.  The  modelling  of  the  decomposition  processes  involves 
sophisticated  numerical  models,  which  require  equations  of  state  describing  properties 
of  the  decomposition  products  such  as  pressure,  volume  and  internal  energy.  These 
decomposition  products  are  largely  gases,  which  in  the  high  pressure,  high 
temperature  regime  present  dining  the  combustion  or  detonation,  are  highly  non-ideal. 
The  non-ideal  behaviour  of  gases  is  largely  determined  by  the  interactions  between 
pairs  of  gas  molecules  which  can  be  described  as  an  intermolecular  potential. 
Consequently,  if  we  can  obtain  accurate  estimates  of  the  intermolecular  potentials 
present,  we  can  accurately  model  the  decomposition  of  the  energetic  materials,  leading 
to  accurate  estimates  of  their  performance  and  decomposition  behaviour. 

Intermolecular  potential  parameters  are  generally  obtained  empirically  by  choosing 
them  to  be  consistent  with  experimental  results.  However,  if  the  intermolecular 
potentials  can  be  obtained  directly  by  computation,  we  can  calculate  their  parameters 
independently  of  experiments,  which  can  then  be  used  to  validate  the  results  of 
experiments  and  also  to  give  a  predictive  capability  for  areas  where  experiments  have 
not,  or  can  not,  be  conducted.  Ab  initio  molecular  orbital  theory  enables  the  calculation 
of  molecular  parameters,  such  as  geometry  and  energy,  using  only  approximations  to 
the  Schrodinger  equation,  without  requiring  any  experimental  data  or  assumptions 
from  the  system  under  study.  The  purpose  of  the  current  study  was  to  attempt  to 
obtain  intermolecular  potentials  for  gaseous  products  of  energetic  material 
decomposition  for  use  in  equations  of  state  using  only  ab  initio  molecular  orbital 
calculations1.  It  was  found  that  these  potentials  can  be  calculated  to  a  good  degree  of 
accuracy,  giving  valuable  inputs  for  equations  of  state  and  allowing  the  subsequent 
calculation  of  the  performance  of  energetic  materials. 


1  The  work  presented  in  this  document  was  commenced  while  Dr  Alex  White  was  on  Long 
Term  Attachment  to  the  Naval  Surface  Warfare  Center,  Indian  Head,  Maryland  and 
subsequently  completed  after  his  return  to  DSTO  Salisbury. 


Authors 


A.  White 

Weapons  Systems  Division 

Alex  White  is  a  Research  Scientist  in  Propulsion  Systems 
Technology  group  of  Weapons  Systems  Division,  DSTO.  He 
obtained  a  Ph.D.  in  1987  in  physical  inorganic  chemistry  studying 
ligand  exchange  processes  in  lanthanide  and  other  complexes.  He 
joined  DSTO  in  1988  and  has  since  worked  in  several  diverse  areas 
including  gun  propellant  service  life,  cost-benefit  analysis  of 
insensitive  munitions,  weapons  effectiveness  studies  and 
computational  chemistry.  He  recently  returned  from  a  long-term 
attachment  to  the  Naval  Surface  Warfare  Center,  Indian  Head, 
MD  in  the  United  States  of  America  where  this  work  was 
commenced. 


F.J.  Zerilli 

Energetic  Materials  Research  and  Technology  Department 
Naval  Surface  Warfare  Center  Indian  Head  Division 
Indian  Head,  MD  20640-5035,  USA 

Frank  J.  Zerilli  is  a  Research  Scientist  in  the  Energetic  Materials 
Research  and  Technology  Department  of  the  Naval  Surface 
Warfare  Center,  Indian  Head,  Maryland,  USA.  He  received  his 
Ph.D.  in  1969  in  theoretical  physics  studying  the  physics  of  black 
holes.  After  several  years  of  teaching  and  research  at  universities, 
he  joined  NSWC  in  1979  where  he  has  done  research  on  the 
equation  of  state  and  strength  of  materials,  both  energetic  and  non- 
energetic. 


H.D.  Jones 

Energetic  Materials  Research  and  Technology  Department 
Naval  Surface  Warfare  Center  Indian  Head  Division 
Indian  Head,  MD  20640-5035,  USA 


Contents 


1.  INTRODUCTION . 1 

2.  THEORY  . 2 

2.1  Energetic  Materials  Decomposition . 2 

2.1.1  Shock  in  an  explosive  without  reaction . 2 

2.1.2  Shock  in  an  explosive  with  reaction  (detonation) . 3 

2.1.3  Combustion  of  propellants . 5 

2.2  Equations  of  State . 5 

2.2.1  The  ideal  gas  law . 6 

2.2.2  The  van  der  Waals  equation  of  state . 6 

2.2.3  The  virial  equation  of  state . 7 

2.2.4  The  Becker-Kistiakowsky-Wilson  (BKW)  equation  of  state . 8 

2.2.5  The  Jones-Zerilli  equation  of  state . 9 

2.3  Intermolecular  Potentials . 10 

3.  TECHNIQUES . 11 

3.1  Introduction . 11 

3.2  Ab  initio  methods . 12 

3.3  Molecular  geometry  optimisation . 13 

3.4  Specifying  molecular  orientation . 14 

3.5  Monte  Carlo  technique . 15 

3.6  Multipole  expansion  technique  for  axially  symmetric  molecules . 16 

4.  RESULTS  AND  DISCUSSION . 18 

4.1  Hydrogen  fluoride  (HF) . 18 

4. 1 .1  Monte  Carlo  technique . 18 

4.1.2  Multipole  expansion  technique . 20 

4.2  Hydrogen  (H2) . 22 

4.2.1  Monte  Carlo  technique . 22 

4.2.2  Multipole  expansion  technique . 23 

4.3  Nitrogen  (N2) . 26 

4.3.1  Monte  Carlo  technique . 26 

4.3.2  Multipole  expansion  technique . 27 

4.4  Carbon  monoxide  (CO) . 29 

4.4.1  Monte  Carlo  technique . 29 

4.4.2  Multipole  expansion  technique . 31 

4.5  Carbon  dioxide  (CO2) . 34 

4.5.1  Monte  Carlo  technique . 34 

4.5.2  Multipole  expansion  technique . 35 

5.  EXAMPLE  -  SHOCK  HUGONIOT  FOR  LIQUID  NITROGEN . 38 

6.  CONCLUSIONS  AND  FURTHER  WORK . 40 

7.  REFERENCES . 40 

APPENDIX  A: ANGULAR  CONFIGURATIONS  FOR  MULTIPOLE  EXPANSION 
TECHNIQUE  . 43 


DSTO-TR-1016 


1.  Introduction 


The  weapons  systems  of  the  Australian  Defence  Force  rely  on  energetic  materials  for 
rocket  and  projectile  propulsion  (propellants)  and  terminal  effects  (explosives).  In 
addition  to  research  aimed  at  decreasing  energetic  material  sensitivity  for  greater 
safety,  improvements  in  performance  are  also  being  sought  for  increased  range  and 
thrust  control  in  propellants  and  for  enhanced  terminal  effects  in  explosives. 

It  is  necessary  to  model  the  decomposition  of  energetic  materials  (combustion  for 
propellants  and  detonation  for  explosives)  in  order  to  gain  a  predictive  capability  for 
the  performance  of  existing  and  potential  energetic  materials  which  may  be  utilised  in 
future  weapons  systems.  This  modelling  involves  sophisticated  numerical  models  of 
the  decomposition  processes  [1],  which  in  turn  require  equations  of  state  describing 
properties  of  the  decomposition  products  such  as  pressure,  volume  and  temperature. 
These  decomposition  products  are  largely  gases,  which  in  the  high  pressure,  high 
temperature  regime  present  during  combustion  or  detonation,  are  highly  non-ideal. 
The  non-ideal  behaviour  of  gases  is  largely  determined  by  the  interactions  between 
pairs  of  gas  molecules  which  can  be  described  by  an  intermolecular  potential. 
Consequently,  if  we  can  obtain  accurate  estimates  of  the  intermolecular  potentials 
present  between  gaseous  products  of  energetic  material  decomposition,  we  can 
realistically  model  the  decomposition  processes  which  occur  in  energetic  materials  to 
obtain  accurate  estimates  of  their  performance  and  decomposition  behaviour. 

Intermolecular  potential  parameters  are  generally  obtained  empirically  by  choosing 
them  to  be  consistent  with  experimental  results.  However,  if  intermolecular  potential 
parameters  can  be  obtained  directly  by  computation  independently  of  experiments, 
they  can  then  be  used  to  validate  the  results  of  experiments  and  also  to  give  a 
predictive  capability  for  areas  where  experiments  have  not,  or  can  not,  be  conducted. 

Ab  initio  molecular  orbital  theory  enables  the  calculation  of  molecular  parameters  such 
as  geometry  and  energy,  using  only  approximations  to  the  Schrodinger  equation, 
without  requiring  any  experimental  data  or  assumptions  from  the  system  under  study. 
The  purpose  of  the  current  study  was  to  attempt  to  obtain  intermolecular  potentials  for 
gaseous  products  of  energetic  material  decomposition  for  use  in  equations  of  state 
using  only  ab  initio  molecular  orbital  calculations.  For  simplicity,  only  axially 
symmetric  (linear)  molecules  will  be  considered  here.  A  more  complex  example  (H2O) 
will  be  addressed  in  a  subsequent  publication  [2]. 

This  work  was  commenced  while  Dr  Alex  White  was  on  Long  Term  Attachment  to  the 
Naval  Surface  Warfare  Center,  Indian  Head  Division,  Indian  Head,  Maryland,  United 
States  of  America. 
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2.  Theory 


2.1  Energetic  Materials  Decomposition 

2.1.1  Shock  in  an  explosive  without  reaction 

A  shock  travelling  through  a  material  such  as  a  solid  explosive  constitutes  a 
discontinuous  (or  approximately  so)  jump  in  the  pressure  and  particle  velocity.  A  solid 
explosive  at  atmospheric  pressure  (Po)  occupies  a  specific  volume  Vo  (Vo  =  1/  po  where 
po  is  the  initial  density  of  the  explosive),  has  specific  internal  energy  Eo  and  initial 
particle  velocity  Uo.  If  a  shock  wave  passes  through  the  material  with  velocity  D  but 
does  not  heat  the  material  sufficiently  to  initiate  reaction,  the  shocked  material  will  be 
at  a  new  specific  volume  V,  pressure  P  and  specific  internal  energy  E.  The  new  particle 
velocity  is  U.  Assuming  steady-state  conditions,  laws  of  conservation  of  energy, 
momentum  and  mass  hold  across  the  shock  (collectively  called  the  "Rankine-Hugoniot 
conditions")  [1,  3, 4, 5]: 


•  conservation  of  energy: 

E-E0=~(P  +  PJ(Ve-V) 


conservation  of  momentum: 


/r-.  tt  \2  t  ,  2  (P  Po)  „2/r\  tt  \2  (P  ^0 ) 

(D~U0)  =V0  - —  or  p0 (D-Uq)  = - — 

V  0/  0  (V0-V)  0  (V0-V) 


conservation  of  mass: 


V  (D-U) 


Vo  (D-U  0) 


(1) 

(2) 


(3) 


Substituting  an  equation  of  state  (see  Section  2.2  below)  describing  E  (or  temperature, 
T)  in  terms  of  V  and  P  into  equation  (1)  gives  a  curve  in  the  V-P  plane  passing  through 
the  initial  state  (Vo,  Po)  called  the  shock  Hugoniot  (Figure  1).  It  is  the  locus  of  all 
possible  states  possible  for  the  explosive  produced  by  shocks  of  varying  strength  from 
that  initial  state  (again  assuming  no  reaction).  Equation  (2)  defines  a  line  plotted  in  the 
V-P  plane  called  the  Rayleigh  line  joining  the  initial  state  (Vo,  Po)  and  the  shocked  state 
(V,  P).  The  slope  of  the  Rayleigh  line  through  (V,  P)  and  (Vo,  Po)  is  -po(D-Uo)2,  and  so  is 
proportional  to  the  square  of  the  shock  velocity  D.  The  Rayleigh  line  represents  the 
jump  condition  from  the  unshocked  to  the  shocked  state  [4]  and,  knowing  (V,  P)  and 
(Vo,  Po),  can  be  used  to  obtain  the  shock  velocity  D  or  alternatively,  knowing  D  and  (Vo, 
Po),  can  be  used  to  calculate  the  shocked  state  (V,  P).  Equation  (3)  can  then  be  used  to 
obtain  U,  the  particle  velocity  behind  the  shock. 
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Figure  1:  Shock  Hugoniot  curve  and  Rayleigh  lines  for  a  solid  explosive  (without  reaction) 

2.1.2  Shock  in  an  explosive  with  reaction  (detonation) 

If  the  shock  in  the  explosive  is  sufficiently  strong,  it  may  heat  the  explosive  sufficiently 
to  cause  a  reaction,  in  which  case  a  reaction  zone  will  form  behind  the  shock  front. 
There  will  be  a  family  of  reaction  Hugoniots,  one  for  each  reaction  coordinate 
extending  through  the  reaction  zone  [5].  The  final  reaction  Hugoniot  is  the  "detonation 
Hugoniot"  or  the  "Hugoniot  of  tide  reaction  products"  (Figure  2).  It  is  displaced  from 
the  shock  Hugoniot  for  the  unreacted  explosive  by  the  heat  of  reaction  (q  =  -AH)  and  so 
does  not  pass  through  the  initial  point  (Vo,  Po).  Behind  the  reaction  zone  the  gaseous 
detonation  products  are  free  to  expand  through  a  rarefaction  ("Taylor")  wave.  The 
products  expand  along  an  isentrope  which  is  closely  approximated  by  the  detonation 
Hugoniot  so  that  it  can  be  considered  that  the  products  "unload"  via  a  family  of 
Rayleigh  lines  in  infinitesimal  steps  down  the  detonation  Hugoniot  [4].  Thus  the  slope 
of  the  detonation  Hugoniot  in  the  shocked  state  (ie,  the  slope  of  the  first  infinitesimal 
unloading  Rayleigh  line)  is  proportional  to  the  square  of  the  velocity  of  the  rear  of  the 
reaction  zone  (or,  equivalently,  the  front  of  the  rarefaction  wave). 

The  conservation  rules  given  above  still  apply,  if  we  assume  that  the  steady-state 
conditions  include  the  reaction  zone,  ie  that  the  rear  of  the  reaction  zone  moves  at  the 
same  speed  as  the  shock  front  and  V  and  P  are  now  the  specific  volume  and  pressure  of 
the  reaction  products  rather  than  the  solid  explosive.  Again,  an  equation  of  state  for  the 
reaction  products  describing  E  (or  T)  in  terms  of  V  and  P  is  required  to  obtain  the 
Hugoniot  for  the  reaction  products. 
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Three  possible  situations  can  be  imagined  regarding  the  relative  position  of  a  particular 

Rayleigh  line  corresponding  to  a  shock  state  in  the  solid  explosive  and  the  detonation 

Hugoniot  as  shown  in  Figure  2: 

(a)  the  Rayleigh  line  may  pass  below  the  detonation  Hugoniot  (as  in  (Vo,  Po)  -  (Vi,Pi)). 
The  shock  velocity  is  insufficient  to  cause  reaction  to  detonation  (ie,  there  cannot  be 
a  steady-state  detonation,  only  a  deflagration). 

(b)  the  Rayleigh  line  may  pass  through  the  detonation  Hugoniot  with  two  points  of 
intersection  (as  in  (Vo,  Po)  -  (V3,P3)).  At  the  upper  intersection  point  (V4,  P4),  the  rear 
of  the  reaction  zone  would  have  a  greater  velocity  (slope  of  the  Hugoniot)  than  the 
shock  front  (slope  of  the  Rayleigh  line)  and  would  overtake  it,  in  violation  of  the 
steady-state  assumption.  At  the  lower  point  (V5,  P5),  the  rear  of  the  reaction  zone 
would  have  a  lower  velocity  than  the  shock  front  and  would  increasingly  lag 
behind  it,  again  violating  the  steady-state  assumption. 

(c)  the  Rayleigh  line  and  the  detonation  Hugoniot  may  touch  tangentially  (as  in  (Vo, 
Po)  -  (V2,P2))  at  the  so-called  Chapman-Jouguet  (CJ)  point.  At  this  point,  the 
velocities  of  the  reaction  zone,  reaction  products  and  the  shock  front  in  the  solid 
explosive  are  all  the  same  (Dq),  and  so  this  is  the  only  point  on  the  detonation 
Hugoniot  where  a  steady-state  detonation  can  occur  [4]. 


Figure  2:  Hugoniot  curves  for  solid  explosive  and  reaction  products  and  Rayleigh  lines 


In  summary,  a  steady-state  shock  front  passing  through  an  explosive  with  velocity  Dq 
compresses  the  solid  explosive  to  a  specific  volume  V2,  raising  the  pressure  to  P2  (and 
raising  the  temperature  to  a  value  determined  by  the  equation  of  state  for  the  solid 
explosive),  thus  initiating  chemical  reaction.  A  reaction  zone  follows  the  shock  front, 
and  on  completion  of  the  reaction,  the  reaction  products  are  at  the  volume  and 
pressure  at  the  CJ  point  on  the  detonation  Hugoniot.  These  products  then  expand  in  a 
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rarefaction  ("Taylor")  wave  behind  the  reaction  zone.  In  a  steady-state  detonation,  the 
shock  front,  reaction  zone  and  front  of  the  rarefaction  wave  are  moving  at  the  same 
velocity  (Dcj). 

Thus,  if  the  elemental  composition,  density  and  heat  of  formation  of  a  solid  explosive 
are  known,  an  equation  of  state  relating  E,  P  and  V  of  the  reaction  products  can  be  used 
in  conjunction  with  the  conservation  relationships  described  above  to  calculate  C-J 
detonation  parameters  of  an  explosive  such  as  equilibrium  detonation  product 
composition  and  pressure  and  temperature.  These  calculations  are  described  in  detail 
in  [1]  using  the  BKW  computer  code. 

2.2.3  Combustion  of  propellants 

Equations  of  state  for  mixtures  of  gases  and  solid  decomposition  products  by  definition 
do  not  depend  on  the  path  taken  to  reach  that  state.  Consequently,  in  addition  to 
detonation  properties  of  explosives,  equation  of  states  are  equally  applicable  to 
prediction  of  propellant  performance. 


2.2  Equations  of  State 

State  functions  are  functions  which  only  depend  on  the  current  state  of  the  system,  not 
on  it's  history.  Examples  include  pressure,  temperature  and  volume.  State  functions  of 
a  particular  system  are  related  by  an  equation  of  state  (EOS).  As  described  above,  an 
EOS  relating  internal  energy  E,  pressure  P  and  volume  V  is  required  to  substitute  into 
the  conservation  equations  given  above  to  solve  them  and  obtain  the  relevant 
Hugoniots  and  detonation  or  combustion  properties. 

The  relationship  of  internal  energy  to  volume  is  given  by  [7]: 


'BE 

[dV 


\ 

Jt 


=  T 


dP) 

dr 


Jv 


-p 


Integrating  this  expression  at  constant  temperature  and  fixed  composition  gives  the 
energy  of  a  system  at  temperature  and  volume  (T,  V)  in  terms  of  its  energy  at 
temperature  T  and  reference  volume  V*  [1]: 


E(T,V)  =  E(T,V*)  + 


dV 


where  (dP I  dT)v  is  obtained  from  the  EOS  relating  P,  V  and  T  (see  below). 
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2.2.1  The  ideal  gas  law 

One  of  the  simplest  and  best  known  equations  of  state  is  the  "ideal  gas  law"  which 
relates  the  pressure  (P),  volume  (V),  temperature  (T)  and  amount  of  gas  (n)  using  the 
universal  gas  constant  (R): 


PV  =  nRT 

or  if  we  define  molar  volume  as  Vm=V/n  : 

PVm  =  RT 

For  an  ideal  gas,  ( dP/dT)v  =nR/V,  so  that  ( dE/dV)T  =  nRT IV  -  P  =  0,  giving 
(dE/dT)p  =CV  =(dE/dT)v  where  Cv  is  the  heat  capacity  of  the  gas  at  constant 
volume.  Thus,  for  an  ideal  gas,  the  internal  energy  is  a  simple  function  of  T,  and  does 
not  depend  on  V  or  P. 

The  ideal  gas  law  assumes  that  there  are  no  forces  present  between  the  molecules  of  the 
gas  and  that  the  gas  molecules  occupy  no  volume.  It  is  only  applicable  where  the  forces 
between  molecules  in  the  gas  are  negligible  compared  with  the  thermal  effects  of  the 
gas,  ie  at  low  pressures.  As  the  pressure  increases,  interactions  between  the  gas 
molecules  become  significant.  At  moderate  pressures,  the  forces  between  the  gas 
molecules  tend  to  be  more  attractive  than  in  the  ideal  gas  case  (where  there  are  no 
forces  between  the  molecules).  This  is  due  to  van  der  Waals  forces  which  arise  from 
electron  correlation  between  molecules1.  At  very  high  pressures,  the  molecules  are 
forced  together  so  closely  that  there  is  non-bonded  overlap  between  the  electron  clouds 
which  in  itself  is  repulsive,  but  also  results  in  incomplete  shielding  of  the  nuclei 
resulting  in  a  strong  electrostatic  repulsion  (short  range).  Consequently,  real  gases  tend 
to  be  easier  to  compress  than  ideal  at  moderate  pressures  due  to  long  range  attractive 
forces  but  harder  to  compress  at  high  pressures  due  to  short  range  repulsive  forces  [7], 
Equations  of  state  for  real  gases  need  to  be  more  complex  than  the  ideal  gas  law. 

2.2.2  The  van  der  Waals  equation  of  state 

Various  improvements  on  the  ideal  gas  law  have  been  proposed  for  real  gases  which 
exhibit  non-ideal  behaviour.  For  example,  if  we  take  account  of: 

(a)  the  fact  that  the  molecules  in  a  gas  take  up  volume,  so  that  when  a  gas  is 
compressed,  the  molecules  only  have  a  volume  of  V-nb  to  move  into  (where  b  is  the 
volume  occupied  by  one  mole  of  molecules)  and 


1  This  can  be  thought  of  in  a  classical  (non-quantum  mechanical)  sense  in  the  following  way: 
Due  to  electron  movement,  the  electric  field  arising  from  the  electrons  in  a  molecule  at  any 
instant  is  not  completely  uniform,  ie.,  there  is  a  net  dipole.  This  instantaneous  dipole  induces  an 
opposite  electric  dipole  in  nearby  molecules  causing  a  net  attraction.  See  [6]  for  more  details. 
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(b)  the  attractive  force  between  molecules  is  proportional  to  the  density  of  molecules 
and  so  reduces  the  observed  pressure  proportionally  to  the  density  (pressure 
reduced  by  an2/V2), 

the  "van  der  Waals"  EOS  can  be  obtained  [7\. 

_  nRT  _  an 2 
~  (V-nb)  V2 

or  substituting  Vm=V/n : 

p_  RT  a 

<Ym-b)  V2 

Although  the  van  der  Waals  EOS  allows  better  approximation  to  real  systems  through 
the  use  of  the  two  coefficients,  a  and  b,  it  fails  at  high  pressures  and  densities  because  it 
assumes  that  the  molecules  have  a  well-defined  diameter  and  does  not  take  into 
account  the  actual  intermolecular  forces  which  lead  to  the  apparent  molecular 
diameter. 

2.2.3  The  virial  equation  of  state 

The  van  der  Waals  EOS  is  a  particular  case  of  the  virial  EOS.  In  general,  the  virial  EOS 
uses  a  power  series  to  model  the  behaviour  of  non-ideal  gases: 

PVm  ,  B(T)  C(T ) 

- 2-  =  l  +  — —  +  ——■  +  ••• 

RT  Vm  V2 

where  B(T)  is  the  second  virial  coefficient,  C(T)  is  the  third,  etc.  It  can  be  seen  that  if  the 
coefficients  are  all  zero,  this  expression  reduces  to  the  ideal  gas  equation.  It  follows 
then  that  the  virial  coefficients  must  arise  from  the  effects  of  the  intermolecular  forces 
in  a  real  gas.  It  can  be  shown  that  B(T)  depends  on  the  intermolecular  potential  (V(r)) 
between  pairs  of  molecules2  [7]: 

B(T )  =  2n Na^  (l  -  e-nr),kT  )r2dr 

where  Na  is  Avogadro's  constant,  k  is  Boltzmann's  constant  and  r  is  the  intermolecular 
distance.  Consequently,  if  we  know  the  form  of  the  intermolecular  pair  potential  V(r) 
for  a  particular  gas,  we  can  evaluate  B(T)  and  obtain  the  virial  EOS  (at  least 
approximately  to  the  second  order)  for  that  gas. 


2  The  third  virial  coefficient  C(T)  depends  on  interactions  between  groups  of  three  molecules 
and  so  on. 
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2.2.4  The  Becker-Kistiakowsky-Wilson  (BKW)  equation  of  state 

The  Becker-Kistiakowsky-Wilson  (BKW)  EOS  is  based  on  a  repulsive  potential  applied 
to  the  virial  EOS  [1].  If  terms  higher  than  third  order  are  neglected,  the  virial  EOS 
becomes: 


PVm  ,  .  B(T )  C(T) 
RT  Vm  Vl 


by  setting  x  =  B(T)/Vm  this  can  be  rewritten  as: 

PV. 


RT 


=  l  +  x+  fix 


or  approximately: 


PVm 

RT 


•  =  1  +  xe 


fa 


If  a  repulsive  potential  of  the  form  U=A/rn  is  used,  then: 


B(T)  = 


K 


K 


(r  +  ©)° 


or  x 


Vm(T  +  &f 


where  a  =  3/n  and  K  A3/n  (the  constant  ©  is  introduced  to  prevent  the  pressure 
tending  to  infinity  as  the  temperature  approaches  zero.  It  is  given  an  arbitrary  value  of 
400  K  [1]).  Thus  the  BKW  EOS  for  a  pure  gas  is  [1]: 


PVm 

_ m 

RT 


l  +  xe^ 


where 


K 

vm(T+ef 


For  a  mixture  of  detonation  products,  the  more  common  form  of  the  BKW  EOS  is  [1]: 

— -  =  l  +  xePx  where  x  =  K — ^ — l—L— 

RT  Vm(T  +  S) 


where  X,  and  h  are  respectively  the  mole  fractions  and  covolumes3  of  the  detonation 
products  and  a,  (3  and  k  are  EOS  constants,  a  is  generally  around  0.5  and  /3  and  K  can 

3  the  geometrical  covolume  is  defined  as  the  volume  occupied  by  a  molecule  rotating  about  its 
centre  of  mass  x  10.46  [1].  This  represents  the  repulsive  potential  of  each  detonation  product  in  a 
similar  way  to  b  in  the  van  der  Waals  equation  of  state.  Although  a  covolume  can  be  calculated 
from  the  molecular  dimensions,  the  best  method  of  evaluating  the  covolume  of  a  product 
molecule  is  from  the  experimental  shock  Hugoniot  [1]. 
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be  obtained  by  iteratively  fitting  experimental  detonation  data.  Generally,  two  sets  of 
parameters  are  used,  one  based  on  RDX  for  most  explosives  and  one  based  on  TNT  for 
explosives  with  large  amounts  of  solid  carbon  in  the  detonation  products  [1]. 

Mader  [1]  describes  in  detail  the  application  of  the  BKW  computer  code  to  predict  the 
performance  of  explosive  mixtures  and  propellants  using  the  BKW  EOS  [8]  for  gaseous 
detonation  products  and  the  Cowan  EOS  [9]  for  solid  detonation  products. 

2.2.5  The  Jones-Zerilli  equation  of  state 

The  EOS  applicable  to  the  work  described  in  this  document  is  detailed  in  [10].  It  is 
based  on  the  formalism  of  Weeks,  Chandler  and  Anderson  [11]  that  repulsive 
intermolecular  forces  dominate  the  characteristics  of  the  fluid  phase.  It  is  assumed  that 
the  intermolecular  interaction  is  spherically  symmetric  and  can  be  expressed  in  the 
form  [10]: 

V(r)  =  u0(r)  +  w(r) 


where  uo(r)  is  the  reference  repulsive  potential: 


fV(r)-V(rc )  r  <rc 
[0  r>rc 


and  w(r)  is  an  attractive  perturbation: 


w(r)  = 


\V(rc) 

jv(r) 


r<rc 
r  >  rc 


where  r  is  the  intermolecular  distance  and  V(r)  is  the  exp-6  potential  with  the 
temperature  dependent  well  depth  as  described  in  Section  2.3  below.  Thus  the 
potential  is  arbitrarily  taken  to  be  divided  at  the  intermolecular  distance  rc4.  Weeks, 
Chandler  and  Andersen  [11]  took  the  breakpoint  rc  to  be  the  intermolecular  distance  at 
the  maximum  well  depth  (r*  below),  however  Ree  [20]  demonstrated  that  by  taking  rc 
to  be  a  variational  parameter,  the  method  could  be  extended  to  the  high-density 
domain  of  interest  in  studies  of  detonation  phenomena. 


4  The  symbol  X  is  generally  used  for  the  breakpoint  of  the  potential,  however  the  symbol  rc  has 
been  used  here  to  avoid  confusion  as  the  symbol  X  is  used  extensively  below  to  refer  to  a 
different  parameter. 
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2.3  Intermolecular  Potentials 


Historically,  the  Lennard-Jones  (LJ)  potential  [12,  13]  has  been  widely  used  to  model 
intermolecular  interactions.  The  LJ  intermolecular  potential  <p  has  a  twelfth  power 
repulsive  term  and  a  sixth  power  attractive  term: 


<t>  =  4e 


K  r  J 


12 


(  , t\6 


\T  ) 


where  r  is  the  intermolecular  separation,  a  is  the  intermolecular  distance  where  the 
potential  is  zero  and  e  is  the  potential  well  depth.  The  form  of  the  LJ  potential  is  shown 
in  Figure  3  along  with  the  parameters  just  described.  The  attractive  tail  arises  from  van 
der  Waals  interaction  due  to  electron  correlations  while  the  strongly  repulsive  core  is 
due  to  nonbonded  overlap  between  electron  clouds  [14, 6]. 


Figure  3:  The  form  of  the  Lennard-Jones  (LJ)  and  the  exponential-6  (exp-6)  potentials. 

However,  the  repulsive  term  in  the  LJ  potential  is  known  to  be  much  too  hard  (ie,  the 
repulsive  core  is  too  steep)  at  small  intermolecular  distances.  A  more  realistic  potential 
which  agrees  with  molecular  beam  scattering  data  and  gives  the  correct  functional 
form  at  small  intermolecular  separations  is  the  exponential-6  (exp-6)  potential  [15]: 
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0  = 


(a -6) 


r 

(  r 

6~ 

6e  r  > -a 

{  r  ) 

where  <(>  is  again  the  potential  at  intermolecular  separation  r,  e  is  again  the  potential 
well  depth,  r*  is  now  die  intermolecular  distance  where  the  potential  is  a  minimum  (ie, 
where  §  =  -e)  and  a  is  a  "steepness  factor".  The  form  of  the  LJ  potential  is  shown  in 
Figure  3  along  with  the  parameters  appropriate  to  this  curve.  It  can  be  seen  that  the 
repulsive  part  of  the  curve  is  indeed  "softer"  than  that  of  the  LJ  potential  (the  steepness 
factor  for  this  curve  has  been  scaled  so  that  the  attractive  terms  of  the  two  curves  are 
identical).  In  this  work,  the  exp-6  potential  was  used  exclusively. 

In  addition,  the  well  depth  e  is  made  temperature  dependent  via  the  expression  [16]: 


e 


f 

V 


N 


to  allow  for  multipole  interactions  which  increase  as  the  temperature  decreases.  This 
effect  should  only  be  important  where  the  molecules  have  strong  multipole 
interactions  and  where  they  do  not  have  dipoles,  X  is  often  assumed  to  be  zero.  The  ab 
initio  techniques  described  here  allow  the  calculation  of  these  parameters  with  no 
assumptions  about  the  system  under  study.  Consequently,  they  afford  a  method  of 
checking  if  this  assumption  is  valid. 


3.  Techniques 


3.1  Introduction 

As  stated  before,  the  object  of  this  work  was  to  use  ab  initio  quantum  chemical  methods 
to  obtain  intermolecular  potential  parameters  for  use  in  equations  of  state  for  the 
products  of  energetic  material  decomposition.  The  intermolecular  potential  is  the 
potential  between  two  molecules  and  as  seen  in  Section  2.3,  varies  with  the 
intermolecular  distance  r.  For  real  molecules  which  are  not  spherically  symmetric,  the 
intermolecular  potential  also  varies  in  a  complex  fashion  with  the  relative  orientation 
of  the  two  molecules.  This  is  generally  overcome  by  specifying  a  "spherically 
symmetric"  potential,  averaged  over  all  possible  relative  molecular  orientations  for 
each  intermolecular  distance. 

As  described  in  Sections  3.5  and  3.6  below,  two  different  techniques  (the  "Monte 
Carlo"  (MC)  and  the  "Multipole  Expansion"  (ME)  techniques)  were  developed  to 
obtain  spherically  symmetric  intermolecular  potentials  and  their  associated 
parameters.  Both  techniques  use  ab  initio  methods  to  obtain  the  intermolecular 
potential  between  two  molecules  with  various  relative  molecular  orientations  at  several 
intermolecular  distances.  The  methods  differ  in  that  while  the  MC  method  averages 
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many  randomly  generated  orientations  at  each  intermolecular  distance,  the  ME 
method  averages  several  specifically  selected  orientations  by  solving  the  set  of 
equations  giving  the  radial  dependence  of  the  potential. 

Consequently,  as  discussed  in  more  detail  below,  the  general  sequence  followed  to 
obtain  the  exp-6  potential  parameters  for  the  intermolecular  potentials  between  two 
molecules  was  as  follows: 

1.  Choose  a  suitable  ab  initio  method  to  be  used  for  the  system  of  interest. 

2.  Obtain  an  optimised  geometry  for  each  molecule  in  isolation,  giving  equilibrium 
internal  bond  lengths  and  angles  and  a  reference  molecular  energy.  The  internal 
geometry  is  fixed  using  these  bond  lengths  and  angles  to  make  the  molecules  rigid. 

3.  For  the  two  molecule  system  with  a  particular  intermolecular  separation,  conduct 
ab  initio  calculations  to  obtain  the  energy  of  the  system  for  several  selected  (ME 
method)  or  many  randomly  generated  (MC  method)  molecular  orientations.  The 
intermolecular  potential  is  this  energy  minus  the  reference  energy  of  the  two 
molecules  in  isolation  (ie,  at  infinite  distance). 

4.  Obtain  the  spherically  symmetric  potential  at  that  intermolecular  distance  by 
averaging  these  intermolecular  potentials  using  the  MC  or  ME  method  as 
appropriate  as  described  below. 

5.  Repeat  steps  3  and  4  for  several  intermolecular  separations  over  a  sufficient  range 
of  the  intermolecular  potential. 

6.  Fit  the  data  using  the  exp-6  potential  curve  to  obtain  the  potential  parameters. 

3.2  Ab  initio  methods 

A  commercial  computer  program  (Gaussian  98  [17])  was  used  to  carry  out  the  ab  initio 
calculations  of  molecular  geometry  and  energies.  Ab  initio  methods  are  discussed  in 
depth  in  [18].  For  the  purpose  of  this  study,  it  is  sufficient  to  note  that  a  particular  ab 
initio  method  is  defined  by  two  components,  a  theoretical  method  and  a  basis  set. 

The  theoretical  method  generally  specifies  the  level  of  electron  correlation  considered. 
In  general,  the  greater  the  level  of  electron  correlation  considered,  the  more  accurate 
the  results,  however  at  a  cost  in  the  use  of  computer  resources  (RAM,  hard  disk  space 
and  computation  time).  The  Moller-Plesset  perturbation  method  [19]  was  used  in  this 
study.  For  the  ME  method,  the  perturbation  is  taken  to  the  fourth  order  (MP4)  whereas 
for  die  MC  method  the  perturbation  was  only  taken  to  the  second  order  (MP2). 
Although  the  latter  is  less  accurate,  it  was  necessary  to  use  it  for  the  large  number  of 
calculations  required  in  the  MC  method. 

The  basis  set  gives  a  mathematical  description  of  the  orbitals  used  to  calculate  the  total 
electronic  wavefunction  of  the  system.  Again,  the  larger  the  basis  set,  the  more  accurate 
the  results,  but  the  greater  the  use  of  computer  resources.  For  the  MC  method,  a 
6-311++G(d,p)  basis  set  was  used.  The  "G"  indicates  that  gaussian  type  orbitals  are 
used.  The  "6-311"  indicates  an  inner  shell  atomic  orbital  composed  of  six  primitive 
gaussian  functions  and  a  valence  shell  composed  of  three  sizes  of  basis  function  for 
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each  valence  orbital.  These  basis  functions  are  composed  of  three,  one  and  one 
primitive  gaussian  functions  respectively.  The  "++"  indicates  the  addition  of  diffuse 
functions  (larger  orbitals)  to  both  hydrogen  and  heavy  atoms,  allowing  for  electron 
density  at  distances  far  from  the  nucleus.  In  the  study  of  intermolecular  potentials, 
these  interactions  are  likely  to  be  important.  The  (d,p)  indicates  that  d-orbital  type 
polarisation  functions  are  added  to  heavy  atoms  while  p-orbital  type  polarisation 
functions  are  added  to  hydrogen  atoms.  These  polarisation  functions  allow  the  orbitals 
to  change  shape  to  allow  for  non-uniform  displacement  of  charge  from  the  nucleus, 
increasing  the  accuracy  of  the  calculation  over  calculations  only  permitted  to  use  the 
functions  available  from  the  filled  orbitals.  For  example,  a  carbon  atom  may  have  d 
polarisation  functions  added,  allowing  a  degree  of  "mixing-in"  of  the  higher  d  orbitals 
to  allow  for  non-uniform  charge  distribution  in  bonded  atoms.  For  the  ME  method, 
which  requires  fewer  calculations  than  the  MC  method,  a  larger  basis  set  could  be 
used,  adding  more  polarisation  functions,  giving  a  6-311++G(3df,  2pd)  basis  set. 

The  computer  resources  required  for  ab  initio  calculations  also  increase  rapidly  with  the 
size  of  the  system  (number  of  atoms)  under  consideration.  Thus  there  is  always  a 
compromise  between  the  accuracy  of  the  result  and  the  computer  resources  and  time 
available. 

3.3  Molecular  geometry  optimisation 

For  both  techniques,  the  first  step  is  to  optimise  the  bond  lengths  and  angles  for  one 
molecule  in  isolation  at  the  level  of  the  subsequent  calculation  (Table  1).  This  is  done  by 
allowing  the  bond  lengths  and  angles  to  change  until  a  minimum  in  the  molecular 
energy  is  found.  The  reference  energy  for  the  two  molecule  system  is  taken  as  twice 
this  energy5. 

Table  1:  Optimised  Molecular  Parameters 


Method 

Theoretical 
Method/Basis  Set 

Species 

Bond 

Length  (A) 

Energy  (hartrees) 

MC 

MP2/6-311++G(d,p) 

HF 

0.916600 

-100.27889234382 

h2 

0.738400 

-1.1603010878136 

n2 

1.120264 

-109.30155631903 

C02 

1.169954 

-188.20655793738 

CO 

1.139923 

-113.07816919596 

ME 

MP4/6-311++G(3df,2pd) 

HF 

0.917878 

-100.34236767 

h2 

0.741210 

-1.1717385902 

n2 

1.113093 

-109.37982637 

C02 

1.172760 

-188.34271091 

CO 

1.142857 

-113.16255765 

5  Alternatively,  the  reference  energy  may  be  taken  as  the  total  energy  of  two  molecules  at 
effectively  infinite  separation.  As  expected,  these  values  are,  for  practical  purposes,  identical. 
However,  this  can  only  be  said  of  methods  which  are  "size-consistent",  ie,  do  not  depend  on  the 
number  of  orbitals  in  fire  system.  MP  methods  are  size-consistent. 
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3.4  Specifying  molecular  orientation 

We  need  to  specify  the  orientation  of  two  molecules  relative  to  one  another  so  that  we 
may  calculate  the  intermolecular  potential  for  that  particular  orientation  and  then 
average  these  potentials  for  several  orientations  using  either  the  Monte  Carlo  method 
(described  in  Section  3.5  below)  or  the  Multipole  Expansion  method  (described  in 
Section  3.6). 

We  assume  that  the  molecules  under  study  are  rigid  bodies,  ie  that  the  bond  lengths 
are  fixed  at  the  bond  length  of  minimum  molecular  energy  for  the  individual  molecules 
obtained  in  Section  3.1.  For  each  molecule,  a  set  of  Cartesian  coordinate  axes  can  be 
assigned  with  the  origin  at  the  centre  of  mass  (CoM)  of  the  molecule.  The  relative 
orientation  of  any  two  molecules  can  then  be  defined  relative  to  an  imaginary  line 
joining  the  centres  of  mass  of  the  two  molecules.  In  general,  this  is  done  by  specifying 
three  angles  ("Euler  angles")  for  each  molecule.  However,  symmetry  considerations 
may  reduce  the  number  of  angles  that  it  is  necessary  to  specify.  For  two  axially 
symmetric  molecules  (ie,  molecules  which  have  a  C«,v  axis),  only  three  angles  need  to 
be  specified  as  shown  in  Figure  4  [22].  The  polar  angles  (3i  and  p2  are  the  angles  that 
molecules  AB  and  CD  make  with  the  intermolecular  axis,  while  ai  is  the  dihedral  angle 
A  -  CoMab  -  CoMcd-  C.  The  molecules  are  separated  by  an  intermolecular  distance  r. 


Figure  4:  Relative  orientation  for  two  axially  symmetric  molecules  AB  and  CD  [22]. 
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3.5  Monte  Carlo  technique 


The  Monte  Carlo  technique  allows  the  calculation  of  temperature  dependent 
intermolecular  potential  parameters.  The  method  for  two  molecules  separated  by  an 
intermolecular  distance  r: 

1.  A  large  number  of  calculations  (1750  in  the  cases  discussed  here)  are  made  over  the 
range  of  possible  relative  molecular  orientations,  values  of  the  internal  angles 
defining  the  relative  intermolecular  orientation  (for  axially  symmetric  molecules, 
the  three  angles  Pi  ,  p2  and  CX2).  In  practice,  this  is  done  by  generating  random 
values  of  Pi  and  p2  from  0  - 180°  and  0.2  from  0  -  360°. 

2.  The  potential  of  the  system  $  is  calculated  for  each  configuration. 

3.  The  configurations  are  averaged  using  summation  as  an  approximation  to  the 
integral: 


0 


( 

-  -kT  In 

v 


1 


J 


4. 

5. 


where  Qi  and  Q2  are  the  orientations  of  the  two  molecules. 

Steps  1  -  3  are  repeated  for  each  value  of  r  over  a  sufficient  range  to  enable  a 
reasonable  curve  fit  of  the  data  around  tire  potential  minimum. 

Temperature  dependent  parameters  e  (maximum  well  depth  of  potential),  r* 
(intermolecular  distance  at  maximum  well  depth)  and  a  (steepness  parameter)  are 
obtained  by  fitting  the  variation  in  potential  with  intermolecular  distance  r  to  the 
exponential-6  (exp-6)  potential: 


(a -6) 


6.  Following  Ree  [20],  the  well  depth  at  infinite  temperature  E0  and  the  temperature 
dependence  of  the  well  depth  X  is  obtained  by  fitting  s  and  T  to: 


7.  Potential  curves  at  different  temperatures  can  be  generated  by  changing  the 
temperature  of  the  calculation.  This  changes  the  average  potential  at  each 
intermolecular  distance  because  of  the  factor  of  kT  in  the  integral. 


The  accuracy  of  the  results  obtained  using  the  Monte  Carlo  method  depends  primarily 
on  three  factors:  the  number  of  points  obtained,  the  theoretical  method  used  and  the 
basis  set  used.  Improvements  in  any  of  these  three  factors  increases  the  computation 
time  required.  It  was  found  that  obtaining  1750  points  at  each  intermolecular  distance 
using  a  second  order  Moller-Plesset  perturbation  method  (MP2)  with  a  6-311++G(d,p) 
basis  set  combination  gave  a  satisfactory  trade-off  between  accuracy  and  computation 
time  for  the  majority  of  the  systems  discussed  here  (see  [18]  for  a  detailed  discussion  of 
ab  initio  methods). 
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3.6  Multipole  expansion  technique  for  axially  symmetric  molecules 

For  two  axially  symmetric  molecules,  approximate  values  for  the  spherical  part  of  the 
potential  may  be  obtained  with  far  fewer  calculations  than  in  the  Monte  Carlo 
technique.  For  these  two  axially  symmetric  molecules,  with  orientations  as  described  in 
Section  3.4,  the  intermolecular  potential  can  be  written: 

V(r, ft, &,«,)  =  jvlo(r,ft,02,a2) 


k=0 


The  functions  yf~k'k)  ,  giving  the  angular  dependence  of  the  potential,  are 

multipole  expansion  functions  of  angles,  in  terms  of  Clebsch-Gordan  coefficients  and 
generalised  spherical  harmonics.  For  example,  the  first  few  terms  in  the  expansion  are: 

Voo=/<Tw 

VJ„  =  cos/5,//1'®  (r)-cos/3,//011  (r) 

V20  =i(3cos2A  +l)//2'0>(r)+i(3cos2ft  +l)/2,0'2,(r) 

+(sin/3]  sin/J2cos a2  -2cosj5j  cos /32)/2CU)(r) 

The  functions  f^l~k'k)  (r)  give  the  radial  dependence  of  the  potential.  The  fta~k’k)(r) 
reduce  at  large  separations  to  the  usual  expansion  in  terms  of  the  electrostatic  spherical 
tensor  multipole  moments  of  molecules  1  and  2: 

f  <*-*.*>  o')  47r  Qi-k,oQk,o 

‘  ^2{l-k)  +  U2k  + 1  r,+1 

For  axially  symmetric  molecules,  an  approximate  result  may  be  obtained  by  solving  the 
set  of  equations  for  )  up  to  fixed  multipole  order.  For  example,  for  1=2 

(quadrupole),  there  are  six  unknown  functions: 

/o(WV),  rw,  /i(W)(r),  /2(2'0)(r),  f2°'2)(r), 


For  1=4  (hexadecapole),  this  increases  to  fifteen: 
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’(3.0)  (  f 


r  ( r),  /3(0'3)  (r),  /4(4'0)  (r),  f™  (r),  /4(2’2)  (r),  /4< 


VX/fV) 


Hence,  if  we  know  V(r,j3l,p2,a2)  for  six  different  angular  configurations  (in  the 
quadrupole  case)  or  fifteen  different  configurations  (in  the  hexadecapole  case)  we  may 
solve  the  set  of  six  (quadrupole)  or  fifteen  (hexadecapole)  linear  equations  for  the 
unknowns  (if  they  are  linearly  independent).  This  number  reduces  if  the  symmetry  of 
the  molecule  renders  some  of  the  unknowns  redundant. 

A  set  of  seventeen  possible  molecular  configurations  is  given  in  the  Appendix,  shown 
for  a  non-identical  pair  of  axially  symmetric  molecules.  Three  linear  combinations  were 
used  to  obtain  the  spherically  symmetric  part  of  the  intermolecular  potential  (V o),  one 
at  the  1=2  (quadrupole)  level  (labelled  CB)  and  two  at  the  1=4  (hexadecapole)  level 
(labelled  CAD  and  CBAD): 

CB:  V0  =  V{C\)  -  —  V  (C2)  +  -V  (C3)  +  -V(Bl)  +  ~  V(B2)  +  \v{B3)  +  \v(B4) 

12  12  2  6  6  6  6 


CAD  :  V0  =  -V(C1)  +  —  7(A1)  +  —  V  (A4)  +  —  V  (Dl)  +  — V(D4) 

0  5  30  30  15  15 

CBAD  :  V0  =  A-V  (Cl)  +  -~V  (C2)  -  A^-V  (C3)  +  j^-V  (Bl)  +  Ty  (B2) 

+  Ay  (S3)  +  —V  (BA)-  —  V  (Al)  -  —V  ( A4)  +  — V  (Dl) 

14  14  105  105  15 

+  —  V(D4)  +  ~V(D2)  +  — V(D5) 

15  15  15 

If  the  molecules  are  identical  (eg,  HF)  several  configurations  become  redundant  (ie,  the 
following  pairs  of  configurations  are  identical:  A4  and  Al,  B4  and  Bl,  B3  and  B2,  D4 

and  Dl,  D5  and  D2,  and  D6  and  D3)  and  the  corresponding  expression  for  the 

spherically  symmetric  part.  If  both  molecules  are  identical  and  both  have  centres  of 
symmetry  (ie,  are  members  of  the  D„h  point  group,  eg,  H2  or  CO2),  the  set  of 
configurations  reduces  even  further  (A1=A2=A3=A4,  B1=B2=B3=B4,  C1=C2,  D4=D1, 
D5=D2,  D6=D3).  Redundancies  in  configurations  for  other  symmetries  (eg,  one  H2 
molecule  and  one  HF)  can  be  analogously  derived. 

For  the  Multipole  Expansion  technique  which  involves  the  gathering  of  far  fewer 
points  compared  to  the  Monte  Carlo  technique,  the  method  and  basis  set  could  be 
improved  in  most  cases  to  MP4/6-311++G(3df,  2pd). 
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4.  Results  and  Discussion 

4.1  Hydrogen  fluoride  (HF) 

4.1.1  Monte  Carlo  technique 

The  average  interxnolecular  potentials  obtained  via  the  Monte  Carlo  (MC)  technique 
and  the  exp-6  curve  fits  for  two  hydrogen  fluoride  (HF)  molecules  at  various 
temperatures  are  shown  in  Figure  5.  The  parameters  obtained  from  the  curve  fitting  are 
given  in  Table  2.  Also  shown  for  comparison  are  results  from  Zerilli  and  Jones  [21] 
derived  from  fitting  molecular  beam  spectroscopic  data  of  Barton  and  Howard  [22], 
labelled  ZJBH.  As  expected,  as  the  temperature  increases,  the  intermolecular  potential 
well  depth  (e)  decreases,  the  radius  at  which  the  well  depth  minimum  occurs  (the 
intermolecular  equilibrium  distance,  r*)  increases  and  the  steepness  parameter  (a)  also 
increases.  The  Monte  Carlo  curve  for  100  000  K  (effectively  infinite  temperature)  is  very 
similar  to  the  ZJBH  curve  at  infinite  temperature,  and  the  associated  parameters 
(Table  2)  are  also  very  similar.  The  corresponding  curve  of  well  depth  variation  with 
temperature  and  the  fit  of  the  relationship  e  =  £o  (1+X/T)  are  shown  in  Figure  6 6  and 
the  parameters  obtained  are  given  in  Table  3.  It  can  be  seen  that  the  Monte  Carlo  and 
the  ZJBH  results  are  generally  in  reasonable  agreement. 


Figure  5:  Average  intermolecular  potential  curves  with  exp-6  fits  (2. 7-3. 6 A)  for  HF/HF  using 
the  Monte  Carlo  technique  (MP2/6-311++G(d,p)) 


18 


6  The  temperature  dependent  well  depth  correction,  e  =  So  (1+A/T),  is  only  valid  at  high 
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Table  2:  HF/HF  Monte  Carlo  (MC)  results  calculated  at  various  temperatures  and  results  of 
Zerilli  and  Jones  (ZJBH)  [21] 


T(K) 

e(K) 

r*  (A) 

a 

MC7 

ZJBH  [21] 

MC3 

ZJBH  [21] 

MC3 

ZJBH  [21] 

300 

1241.9 

1389 

2.80 

2.62 

10.96 

11.58 

2  000 

313.3 

389 

3.04 

3.00 

12.92 

12.24 

3  000 

249.5 

317 8 

3.10 

3.06 

13.11 

12.53 

10  000 

167.0 

202 

3.20 

3.20 

13.43 

13.24 

100  000  (inf) 

138.4 

1529 

3.26 

3.284 

13.60 

13.654 

T  (K) 


Figure  6:  Well  depth  variation  with  temperature  for  HF/HF  using  the  Monte  Carlo  technique 
and  fit  to  £=£o  (1+Z/T)  (fit  700  - 100  000  K) 


Table  3:  Well  depth  relationship  parameters  for  HF/HF 


Source 

eo  (K)  X  (K) 

MC10 

ZJBH  [21] 

169  2090 

152  3250 

7  exp-6  fit  of  2.7  -  3.6  A  Monte  Carlo  data. 

8  from  e=e0(l+X/T). 

9  calculated  at  infinite  temperature. 

10  e=e0(l+X/T)  fit  of  700  - 100  000  K. 
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4.1.2  Multipole  expansion  technique 

The  intermolecular  potential  curves  obtained  for  several  configurations  of  the  HF/HF 
system  are  shown  in  Figure  7. 


r(A) 

Figure  7:  Intermolecular  potential  curves  for  HF/HF  for  several  configurations  used  in  the 
multipole  expansion  technique  (MP4/6-311++G(3df2pd)) 


It  can  be  seen  that  the  potentials  vary  as  might  be  expected  by  examination  of  the 
configurations  given  in  the  Appendix.  For  example,  the  A1  configuration  is  one  in 
which  the  hydrogen  fluoride  molecules  are  arranged  end-on  with  the  hydrogen  of  one 
molecule  facing  the  fluorine  of  the  other  (ie,  H-F—H-F).  Thus,  the  dipoles  of  the 
molecules  are  aligned  and  as  the  molecules  are  brought  towards  one  another  (reducing 
r),  there  is  a  deep  potential  well  at  the  equilibrium  separation  distance  of  about  2.8  A. 
The  potential  rises  rapidly  as  r  is  reduced  further  and  repulsive  forces  predominate.  On 
the  other  hand,  in  the  B1  configuration,  the  hydrogen  fluoride  molecules  are  arranged 
side-on  with  the  dipoles  of  each  molecule  aligned: 

H - F 

I  l 

l  l 

I  I 

I  I 

H - F 

As  the  molecules  are  brought  closer  together,  repulsive  forces  predominate  throughout 
and  there  is  no  minimum  to  the  potential  energy  curve  (or,  more  correctly,  the 
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minirrmm  is  zero  at  r  =  «>).  The  potential  curves  for  the  other  configurations  can  be 
rationalised  in  a  similar  manner. 

The  three  average  potentials  from  the  multipole  expansion  technique  described  in 
Section  3.6  (labelled  ME(CB)/  ME(CAD)  and  ME(CBAD)),  are  shown  in  Figure  8  for  the 
system  of  two  hydrogen  fluoride  molecules.  For  comparison,  also  shown  in  Figure  8 
are  the  curve  (at  effectively  infinite  temperature,  ie  100  000  K)  from  the  Monte  Carlo 
analysis  above  (labelled  MC(inf))  and  the  curve  from  the  data  of  Zerilli  and  Jones  [21] 
calculated  at  infinite  temperature  determined  from  data  of  Barton  and  Howard  [22] 
(labelled  ZJBH(inf)).  The  derived  parameters  are  given  in  Table  4.  It  can  be  seen  that 
the  averages  at  the  hexadecapole  level  (ME(CAB)  and  ME(CBAD))  are  much  more 
accurate  than  that  at  the  quadrupole  level  (ME(CB)).  These  higher  level  averages  are 
close  to  the  ZJBH(inf)  values,  however,  in  this  case,  the  Monte  Carlo  results  are  much 
closer. 


Figure  8:  Average  intermolecular  potential  curves  for  HF/HF  using  the  multipole  expansion 
technique  (MP4/6-311++G(3df,2pd)) 


Table  4:  Summary  of  intermolecular  potential  parameters  for  HF/HF 


Source 

eo  (K) 

r*  (A) 

a 

ME  (CB)11 

504.0 

2.94 

10.8 

ME  (CAD)7 

172.2 

3.19 

12.3 

ME  (CBAD)7 

193.1 

3.16 

12.1 

MC  (inf)12 

138.4 

3.26 

13.6 

ZJBH  (inf)  [21] 

152.0 

3.28 

13.7 

11  exp-6  fit  of  2.5  -  5.0  A  multipole  expansion  data. 

12  exp-6  fit  of  2J  -  3.6  A  Monte  Carlo  data  at  100  000  K  (ie,  effectively  infinite  temperature). 
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4.2  Hydrogen  (H2) 

4.2.1  Monte  Carlo  technique 

Plots  of  the  intermolecular  potential  at  various  temperatures  using  the  Monte  Carlo 
technique  for  the  H2/H2  system  are  shown  in  Figure  9.  Note  the  scale  used  on  the 
ordinate  axis  of  this  graph.  It  can  be  seen  that  there  is  negligible  difference  between  the 
potential  curves  as  the  temperature  increases.  For  comparison,  also  shown  are  results 
from  liquid-state  perturbation  theory  calculations  matched  to  high-pressure  data 
(labelled  JZ)  [10, 24].  The  calculated  parameters  are  given  in  Table  5. 


Figure  9:  Average  intermolecular  potential  curves  with  exp-6  fits  (2.7 -4. 2 A)  for  H2/H2  using  the 
Monte  Carlo  technique  (MP2/6-311++G(d,p)) 


Table  5:  H2/H2  Monte  Carlo  (exp-6  fit  of  2.7  -  4.2  A)  results  calculated  at  various  temperatures 
and  results  of  Jones  and  Zerilli  (JZ)  [10, 24] 


T  (K) 

e(K) 

r*(A) 

a 

300 

10.9 

3.92 

12.9 

2  000 

10.8 

3.92 

12.9 

3  000 

10.8 

3.92 

12.9 

10  000 

10.8 

3.92 

12.9 

100  000  (inf) 

10.8 

3.92 

12.9 

JZ  [10,24] 

36.0 

3.46 

11.1 

The  well  depth  variation  with  temperature  is  given  in  Figure  10.  There  is  negligible 
change  in  the  well  depth  over  the  complete  temperature  range  (noting  again  the  scale 
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on  the  ordinate  axis  and  remembering  that  the  ordinate  scale  in  the  corresponding 
figure  for  the  HF/HF  system  above  (Fig.  6)  was  logarithmic).  The  parameters  obtained 
for  the  H2/H2  system  are  given  in  Table  6.  The  very  small  temperature  dependence 
observed  in  the  graph  is  also  reflected  in  the  very  small  value  of  X. 


Figure  10:  Well  depth  variation  with  temperature  for  H2JH2  using  the  Monte  Carlo  technique 
and  fit  to  £=£o  (1+XfT)  (fit  700  -  100  000  K) 

Table  6:  Well  depth  relationship  parameters  for  H2/H2 


Source 

£0  (K)  X  (K) 

MC13 

JZ  [10, 24] 

10  1 

36  0 

4.2.2  Multipole  expansion  technique 

The  intermolecular  potential  curves  for  H2/H2  for  some  of  the  configurations  from  the 
multipole  expansion  (ME)  technique  are  shown  in  Figure  11.  The  same  curves  are 
shown  in  Figure  12  with  the  ordinate  scale  expanded  to  better  illustrate  the  angular 
dependence  of  the  potential.  It  can  be  seen  that  the  angular  dependence  is  much  less 
significant  for  the  H2/H2  system  than  previously  seen  for  HF/HF  (Fig.  7).  This  is,  of 
course,  to  be  expected  from  the  fact  that  HF  has  a  strong  dipole  moment  which  is 
absent  in  H2. 


13  e=e0(l+^T)  fit  of  700  - 100  000  K. 
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Figure  11:  Intermolecular  potential  curves  for  H2/H2  for  several  configurations  used  in  the 
multipole  expansion  technique  (MP4/6-311++G(3df,2pd)) 


Figure  12:  Intermolecular  potential  curves  (expanded  ordinate  scale)  for  Fh/Hifor  several 

configurations  used  in  the  multipole  expansion  technique  MP4/6-311++G(3df,2pd)) 
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The  plots  of  average  potential  are  shown  in  Figure  13  and  the  parameters  obtained  are 
given  in  Table  7.  Again,  the  corresponding  results  from  the  Monte  Carlo  analysis  at 
100  000  K  (labelled  MC  (inf))  and  results  from  the  data  of  Jones  and  Zerilli  [10,  24] 
(labelled  JZ)  are  given  for  comparison. 


Figure  13:  Average  intermolecular  potential  curves  for  H2/H2  using  the  multipole  expansion 
technique  (MP4/6-31 1  ++G( 3df,2pd)) 


Table  7:  Summary  of  intermolecular  potential  parameters  for  H2/H2 


Source 

80  (K) 

r  (A) 

a 

ME  (CB)14 

29.3 

3.55 

11.8 

ME  (CAD)10 

21.9 

3.65 

12.1 

ME  (CBAD)10 

22.5 

3.64 

12.1 

MC  (inf)15 

10.8 

3.92 

12.9 

JZ  [10,  24] 

36.0 

3.46 

11.1 

It  can  be  seen  that  the  results  at  the  hexadecapole  level  (CAD  and  CBAD)  are  very  close 
to  one  another.  In  this  case,  the  Multipole  Expansion  results  are  closer  than  the  Monte 
Carlo  results  to  those  of  Jones  and  Zerilli  [10, 24]  (in  fact,  the  results  at  the  quadrupole 
level,  CB,  appear  to  be  closest).  However,  given  the  very  small  values  of  £0,  all  the 
results  are  reasonably  close. 


14  fit  2.5 -5.0  A. 

15  exp-6  fit  of  2.7  -  4.2  A  Monte  Carlo  data  at  100  000  K  (ie,  effectively  infinite  temperature). 
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4.3  Nitrogen  (N2) 

4.3.2  Monte  Carlo  technique 

Plots  of  the  intermolecular  potential  for  the  N2/N2  system  at  various  temperatures  are 
shown  in  Figure  14  along  with  the  curve  from  the  data  of  Jones  and  Zerilli  [10, 24].  The 
corresponding  parameters  are  given  in  Table  8.  It  can  be  seen  that  the  attractive  parts  of 
the  potential  curves  are  all  similar  (ie,  not  dependent  on  temperature)  and  that  the 
minima  (well  depths)  of  the  potential  curves  are  all  similar.  This  is  also  reflected  in  the 
plot  of  well  depth  against  temperature  (Fig.  15)  and  in  the  calculated  well  depth 
relationship  parameters  given  in  Table  9. 


Figure  14:  Average  intermolecular  potential  curves  with  exp-6  fits  (3.4-47 A)  for  Nfflzusing 
the  Monte  Carlo  technique  (MP2/6-3 1 1  ++G(d,p)) 


Table  8:  N2/N2  Monte  Carlo  (exp-6  fit  of  3.4  -  4.7  A)  results  calculated  at  various  temperatures 
and  results  of  Jones  and  Zerilli  (JZ)  [10,  24] 


T  (K) 

e(K) 

r*(A) 

a 

300 

223.7 

4.01 

9.9 

2  000 

217.7 

4.08 

11.6 

3  000 

217.9 

4.08 

11.9 

10  000 

218.8 

4.08 

12.6 

100  000  (inf) 

219.4 

4.08 

12.9 

JZ  [10,  24] 

100.0 

4.12 

13.3 
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Figure  15:  Well  depth  variation  with  temperature  for  N2/N2  using  the  Monte  Carlo  technique 
and  fit  to  e  =  so  (1 +X/T)  (fit  700  - 100  000  K) 


Table  9:  Well  depth  relationship  parameters  for  N2/N2 


Source 

eo  (K)  X  (K) 

MO* 

JZ  [10] 

220  27 

100  0 

43.2  Multipole  expansion  technique 

The  intermolecular  potential  curves  for  each  configuration  are  shown  in  Figure  16  and, 
with  an  expanded  ordinate  scale,  in  Figure  17.  The  observed  potential  curves  show  a 
great  deal  of  angular  dependence,  with  the  end-on  configuration  (Al)  rising  rapidly  as 
the  nitrogen  atoms  are  brought  closer  together.  The  average  potential  curves  are  shown 
in  Figure  18  and  the  parameters  obtained  from  the  curve  fitting  are  given  in  Table  10. 


1(5  e  =  e0  (1+A./T)  fit  of  700  - 100  000  K. 
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Figure  16:  Intermolecular  potential  curves  for  N2/N2  for  several  configurations  used  in  the 
multipole  expansion  technique  (MP4/6-311++G(3df2pd)) 


Figure  1 7:  Intermolecular  potential  curves  (expanded  ordinate  scale)  for  N1/N2  for  several 
configurations  used  in  the  multipole  expansion  technique  (MP4/6-311++G(3df2pd)). 
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Figure  18:  Average  intermolecular  potential  curves  for  N2/N2  using  the  multipole  expansion 
technique  (MV 4/6-31 1++G(3df,2pd)). 


Table  10:  Summary  of  intermolecular  potential  parameters  for  N2/N2 


Source 

eo  (K) 

”  ^A) 

a 

ME  (CB)17 

205.4 

4.00 

13.7 

ME  (CAD)i3 

166.2 

4.17 

13.6 

ME  (CBAD)i8 

238.2 

3.85 

13.2 

MC  (inf)i8 

219.4 

4.08 

12.9 

JZ  [10, 24] 

100.0 

4.12 

13.3 

4.4  Carbon  monoxide  (CO) 

4.4.2  Monte  Carlo  technique 

Plots  of  the  intermolecular  potential  for  the  CO/ CO  system  at  various  temperatures  are 
shown  in  Figure  19  along  with  the  curve  from  tire  data  of  Jones  and  Zerilli  [10,  24],  The 
corresponding  parameters  are  given  in  Table  11.  It  can  be  seen  that  the  attractive  parts 
of  the  potential  curves  at  large  intermolecular  distances  and  the  minima  (well  depths) 
of  the  potential  curves  are  all  similar  (ie,  not  dependent  on  temperature).  This  is  also 
reflected  in  the  plot  of  well  depth  against  temperature  (Fig.  20)  and  in  the  calculated 
well  depth  relationship  parameters  given  in  Table  12. 


17  fit  3.0 -5.0  A. 

18  exp-6  fit  of  3.4  -  4.7  A  Monte  Carlo  data  at  100  000  K  (ie,  effectively  infinite  temperature). 
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Figure  19:  Average  intermolecular  potential  curves  with  exp-6  fits  (3. 2-4.6 A)  for  CO/CO  using 
the  Monte  Carlo  technique  (MP2/6-311++G(d,p)) 


Table  11:  CO/CO  Monte  Carlo  (exp-6  fit  of  3.2  -  4.6  A)  results  calculated  at 
temperatures  and  results  of  Jones  and  Zerilli  (JZ)  [10, 24] 


various 


T  (K) 

e(K) 

r*(A) 

a 

300 

129.5 

4.19 

10.0 

2  000 

120.6 

4.39 

11.0 

3  000 

120.0 

4.39 

11.4 

10  000 

118.2 

4.39 

12.4 

100  000  (inf) 

116.4 

4.37 

13.3 

JZ  [10,  24] 

100.0 

4.12 

13.3 
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Figure  20:  Well  depth  variation  with  temperature  for  CO/CO  using  the  Monte  Carlo  technique 
and  fit  to  s=  So  (1+2./T)  (fit  700  -  100  000  K) 

Table  12:  Well  depth  relationship  parameters  for  CO/CO 


Source 

eo  (K)  A  (K) 

MC19 

JZ  [10] 

116  34 

100  0 

4.4.2  Multipole  expansion  technique 

The  intermolecular  potential  curves  for  each  configuration  are  shown  in  Figure  21,  and 
with  an  increased  ordinate  scale  in  Figure  22.  As  in  the  case  of  N2,  there  is  a  great  deal 
of  angular  dependence  of  the  potential  however,  whereas  in  the  N2/N2  case  (Figs.  16 
and  17)  the  curve  for  the  B1  configuration  rises  rapidly  in  a  similar  manner  to  the  curve 
for  the  D1  configuration,  in  the  CO/CO  case,  the  B1  curve  rises  much  slower,  ie,  the 
repulsion  is  lower  at  similar  intermolecular  distances.  This  is  probably  due  to  the 
dipole  present  in  CO,  which  gives  a  8+  charge  to  the  carbon  and  a  8-  charge  to  the 
oxygen  in  each  molecule.  Then,  due  to  the  different  atomic  masses  of  the  carbon  and 
oxygen  atoms,  the  centre  of  mass  of  each  molecule  is  closer  to  the  oxygen  than  the 
carbon.  Consequently,  in  the  B1  configuration  (Fig.  23),  the  carbon  of  one  CO  molecule 
is  close  to  the  oxygen  of  the  other  molecule,  giving  rise  to  an  electrostatic  attraction 
which  counteracts  the  repulsion  from  the  non-bonding  overlap  of  the  electron  clouds. 


19  e  =  e0  (1+A/T)  fit  of  700  - 100  000  K. 
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21:  Intermolecular  potential  curves  for  CO/CO  for  several  configurations  used  in  the 
multipole  expansion  technique  ( MP4/6-31 1  ++G(  3df2pd)) 


Figure  22:  Intermolecular  potential  curves  (expanded  ordinate  scale)  for  CO/CO  for  several 
configurations  used  in  the  multipole  expansion  technique  (MP4/6-311++G(3df2pd)). 
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Figure  23:  B1  configuration  for  CO /CO. 

The  average  potential  curves  are  shown  in  Figure  24  and  the  parameters  obtained  from 
the  curve  fitting  are  given  in  Table  13.  As  might  be  expected,  the  results  are  very 
similar  to  those  for  the  N2/N2  system. 


Figure  24:  Average  intermolecular  potential  curves  for  CO/CO  using  the  multipole  expansion 
technique  (MP4/6-311++G(3df,2pd)). 


Table  13:  Summary  of  intermolecular  potential  parameters  for  CO/CO 


Source 

£0  (K)  r  (A)  a 

ME  (CB)2o 

ME  (CAD)16 

ME  (CBAD)16 

197.1  4.14  13.0 

154.7  4.27  13.3 

205.0  3.99  12.7 

MC  (inf)21 

JZ  [10, 24] 

116.4  4.37  13.3 

100.0  4.12  13.3 

2°  fit  3.0  -5.0  A. 

21  exp-6  fit  of  3.2  -  4.6  A  Monte  Carlo  data  at  100  000  K  (ie,  effectively  infinite  temperature). 


33 


DSTO-TR-1016 


4.5  Carbon  dioxide  (CO2) 

4.5.1  Monte  Carlo  technique 

Plots  of  the  intermolecular  potential  for  the  CO2/CO2  system  at  various  temperatures 
are  shown  in  Figure  25  along  with  the  curve  from  the  data  of  Jones  and  Zerilli  [10,  24]. 
The  corresponding  parameters  are  given  in  Table  14.  Again,  it  can  be  seen  that  the 
attractive  parts  of  the  potential  curves  are  all  similar  (ie,  not  dependent  on 
temperature)  and  that  the  minima  (well  depths)  of  the  potential  curves  are  all  similar. 
This  is  also  reflected  in  the  plot  of  well  depth  against  temperature  (Fig.  26)  and  in  the 
calculated  well  depth  relationship  parameters  given  in  Table  15. 


Figure  25:  Average  intermolecular  potential  curves  with  exp-6  fits  (3.4  -5.0  A)  for  CO2/CO2 
using  the  Monte  Carlo  technique  (MP2/6-311++G(d,p)) 


Table  14:  CO2/CO2  Monte  Carlo  (exp-6  fit  of  3.4  -  5.0  A)  results  calculated  at  various 
temperatures  and  results  of  Jones  and  Zerilli  (JZ)  [10,  24] 


T  (K) 

e(K) 

r'(A) 

a 

300 

329.0 

4.07 

9.4 

2  000 

210.9 

4.57 

10.4 

3  000 

206.4 

4.66 

10.6 

10  000 

216.4 

4.80 

11.5 

100  000  (inf) 

212.5 

4.73 

15.1 

JZ  [10, 24] 

250.0 

4.15 

13.4 

34 


DSTO-TR-1016 


Figure  26:  Well  depth  variation  with  temperature  for  CO2JCO2  using  the  Monte  Carlo 
technique  and  fit  to  e=so  (1+XfT)  (fit  700  - 100  000  K) 


Table  15:  Well  depth  relationship  parameters  for  CO2/CO2 


Source 

eo  (K)  X  (K) 

MC22 

JZ  [10] 

210  144 

250  0 

4.5.2  Multipole  expansion  technique 

Attempts  were  made  to  obtain  EOS  parameters  for  CO2  using  the  multipole  expansion 
technique.  Studies  at  the  MP4/6-311++G(3df,2pd)  level  proved  intractable  on  our 
computers  so  the  studies  were  conducted  at  the  MP2/6-311++G(3df,2pd)  level.  The 
potential  curves  obtained  from  these  studies  for  each  configuration  are  shown  in 
Figure  27.  The  same  data  shown  in  a  plot  with  an  expanded  ordinate  scale  are  shown 
in  Figure  28.  It  can  be  seen  that  there  is  a  significant  amount  of  bonding  interaction 
occurring  at  intermolecular  distances  of  around  3.5  A  in  the  A1  (end-on)  configuration. 
This  is  not  surprising,  as  for  two  CO2  molecules,  when  the  intermolecular  distance 
(distance  between  the  carbon  atoms  on  the  two  different  molecules)  is  3.5A,  the 
distance  between  the  two  closest  oxygen  atoms  is  comparable  to  the  C-O  bond  length 
(1.167 A).  Consequently,  expressions  for  the  spherically  symmetric  potential  containing 
contributions  from  the  A  configurations  (ie,  CAD  and  CBAD)  are  probably  unreliable. 


22  8  =  so  (1+X/T)  fit  of  700  - 100  000  K. 
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Figure  27:  Intermolecular  potential  curves  for  CO2/CO2  for  several  configurations  used  in  the 
multipole  expansion  technique  (MP2/6-311++G(3df,2pd)).  A  bonding  interaction  is 
apparent  in  the  curve  for  the  A1  configuration  at  about  3.5  A. 


Figure  28:  Intermolecular  potential  curves  (increased  ordinate  scale)  for  CO2/CO2  for  several 
configurations  used  in  the  multipole  expansion  technique 
(MP2/6-311++G(3df,2pd)). 
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The  curves  from  the  three  averaging  expressions  are  shown  in  Figure  29  and  the  results 
are  given  in  Table  16.  Again,  the  results  for  CAD  are  probably  unreliable  and  are  given 
for  information  only.  A  meaningful  curve  fit  could  not  even  be  obtained  for 
ME  (CBAD).  The  results  for  CB  should  be  reliable,  however  as  seen  before,  results  from 
this  second  order  multipole  expansion  are  generally  not  as  accurate  as  those  using  a 
fourth  order  expansion. 


Figure  29:  Average  intermolecular  potential  curves  for  CO2/CO2  using  the  multipole 
expansion  technique  (MP2/6-311++G(3df,2pd)).  ME  (CAD)  and  ME  (CBAD) 
data  are  probably  unreliable  due  to  bonding  interactions  in  the  A1  (end-on) 
configuration 


Table  16:  Summary  of  intermolecular  potential  parameters  for  CO7/CO2 


Source 

so  (K) 

r(A) 

a 

ME  (CB)23 

519.9 

4.10 

15.3 

ME  (CAD)24 

139.8 

5.02 

18.4 

ME  (CBAD) 

n/a 

n/a 

n/a 

MC  (inf)25 

212.5 

4.73 

15.1 

TZ  riO,  241 

250.0 

4.15 

13.4 

23  fit  3.5 -6.0  A. 

24  fit  3.8  -  6.0  A.  ME  (CAD)  results  probably  unreliable  due  to  bonding  interactions  as  discussed 
in  text. 

25  exp-6  fit  of  3.4  -  5.0  A  Monte  Carlo  data  at  100  000  K  (ie,  effectively  infinite  temperature). 
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5.  Example  -  Shock  Hugoniot  for  Liquid  Nitrogen 


The  shock  Hugoniot  for  liquid  nitrogen  was  calculated  using  the  parameters  from  the 
multipole  expansion  technique  given  in  Table  10.  The  results  compared  with 
experimental  data  given  by  Dick  [25]  are  shown  in  Figures  30  and  31.  It  can  be  seen 
that  the  calculated  curves  are  a  very  good  match  to  the  experimental  data. 


N2  Hugoniot,  liquid,  T0  =  75  K,  V0  =  1 .21 95  cc/gm,  P0  =  1  atm. 


Figure  30:  Shock  Hugoniot  for  liquid  nitrogen  (pressure  vs  specific  volume) 
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PARTICLE  VELOCITY  (km/s) 


Figure  31:  Shock  velocity  vs  specific  volume  for  liquid  nitrogen 
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6.  Conclusions  and  Further  Work 

This  document  describes  two  methods  for  ab  initio  calculation  of  intermolecular 
potential  parameters  for  use  in  equations  of  state  based  on  the  formalism  of  Week, 
Chandler  and  Andersen  [11]  which  describes  fluids  whose  molecules  interact  via  a 
spherically  symmetric,  modified  Buckingham  (exp-6)  potential.  The  two  methods 
described  are: 

•  a  Multipole  Expansion  method,  suitable  for  analytical  calculation  of  intermolecular 
potential  parameters  for  axially  symmetric  molecules 

•  a  Monte  Carlo  method,  which  can  be  used  to  obtain  average  potential  energy 
parameters  for  any  molecules  within  the  accuracy  of  the  ab  initio  method  used. 

The  results  obtained  are  in  reasonable  agreement  with  those  obtained  from  fitting 
experimental  results  [10,  24].  It  remains  now  to  use  these  parameters  in  calculations  of 
explosive  performance  parameters.  Preliminary  results  for  the  fluorine-containing 
explosive  FEFO  (bis(2-fluoro-2,2-dinitroethyl)formal)  appear  promising  [26]. 

The  systems  under  study  here  constitute  intermolecular  potentials  for  pairs  of  axially 
symmetric  identical  molecules.  For  mixed  molecules  (eg,  HF/H2),  the  Lorentz- 
Berthelot  mixing  rule  is  often  invoked  to  estimate  the  intermolecular  potential  from  the 
potentials  of  the  pairs  of  identical  molecules.  The  methods  described  here  afford  a 
method  of  directly  measuring  these  heteromolecular  potentials  and,  consequently,  of 
checking  the  validity  of  the  parameters  used  in  the  Lorentz-Berthelot  rule.  This  aspect 
will  be  addressed  in  a  forthcoming  publication  [6].  The  application  of  the  Monte  Carlo 
method  to  non-axially  symmetric  molecules,  in  particular  H2O,  will  also  be  described 
in  a  subsequent  publication  [2]. 
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Appendix  A:  Angular  Configurations  for  Multipole 

Expansion  Technique 


D1  (45,45,0)  D2  (45,45,180)  D3  (45,45,90)  D4  (135,135,0)  D5  (135,135,180)  D6  (135,135,90) 
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